General theory

• Estimation: minimizes cost (or maximizes similarity)
• Similarity metrics/cost functions

  • Least squares (images have the same scale) C\{ R(v), I(T(v))\} = \int_{v \subset \mathbb{R}^3} \{ I(T(v)) - R(v)\}^2 dv

  • Correlation (images on different scales)

  • Mutual information (images with different contrasts/distributions) \mathrm{MI}(R, I) = \sum_{r,i} \mathbb{P}_{RI}(r, i) \log \left( \frac{\mathbb{P}_{RI}(r, i)}{\mathbb{P}_R(r), \mathbb{P}_I(i)} \right)

    where \mathbb{P}_{RI}(r, i): is the joint probability, and \mathbb{P}_R(r), \mathbb{P}_I(i) are the marginal probabilities

• Transformation computed to minimize \min_T C\{R(v), I(T(v))\}
• Some of these cost functions are convex given some assumptions on the images
• That would imply a global optimum, given some assumtions on I and R
• Those assumptions are not met, so many adhoc approaches are implemented¹
• Penalty terms are often useful
• Registration is a challenging problem

Rigid-body

• Six parameter rigid registration
• Used to align functional images over time (motion correction)
• Also used to align images from the same participant collect in the same scan session.
• y = (y_1, y_2, y_3, 1) - template coordinates, x = (x_1, x_2, x_3, 1) input coordinates
• Appended with an intercept for transformation
• Can be described as the composition of four transformations

\begin{bmatrix} y_1 \\ y_2 \\ y_3\\ 1 \end{bmatrix} = T_4 T_3 T_2 T_1 \cdot \begin{bmatrix} x_1 \\ x_2 \\ x_3\\ 1 \end{bmatrix}

Translation of image T_4 = \begin{bmatrix} 1 & 0 & 0 & q_1 \\ 0 & 1 & 0 & q_2 \\ 0 & 0 & 1 & q_3\\ 0 & 0 & 0 & 1 \end{bmatrix}

Rotation around the x-axis:

T_1 = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & \cos(q_4) & \sin(q_4) & 0 \\ 0 & -\sin(q_4) & \cos(q_4) & 0\\ 0 & 0 & 0 & 1 \end{bmatrix}

Other rotations are defined similarly. See the reference above.

Interpolation

• Images are measured on a grid, so applying the transformation T may require evaluation of the image, I, at locations not represented by a grid coordinate
• Interpolation computes I(T(v)) at non-grid locations
• Required for evaluating cost function and when applying the resulting transformation matrix.
• Interpolation methods
• Nearest neighbor interpolation is important when registering atlas images
  • Atlas images are integer valued and nearest neighbor ensures output is integer valued

Nonlinear (deformable) registration

• Deformable registration allows different parts of the image to have a unique transformation matrix
• Called a deformation field \phi(x) such that y = \phi(x).
• Different approaches:
  • B‑splines, diffeomorphisms (e.g., LDDMM, SyN), demons, finite-element
• Enforces smoothness, invertibility, and topology preservation via regularization.
• Typically linear transformation is estimated first, then function composition is used when estimating and evaluating the nonlinear transformation
• Composition avoids repeated interpolation of the image
• ANTs registration paper Avants et al (2009).

fnirt --ref=images/MNI152_T1_1mm_brain.nii.gz --in=images/simonReg/simon_brain.nii.gz --aff=images/simonReg/affineReg.txt --fout=images/simonReg/field.nii.gz --jout=images/simonReg/jacobian.nii.gz --iout=images/simonReg/simon_nonlinearReg.nii.gz --cout=images/simonReg/simon2MNI.nii.gz 

Visualize results

fsleyes images/MNI152_T1_1mm_brain.nii.gz images/simonReg/simon_brain_affineReg.nii.gz images/simonReg/simon_nonlinearReg.nii.gz

Derivatives of the registration process

• Features of the registration process are informative of brain structure
• Probabilistic tissue segmentation: images can be analyzed across participants after registration to template space
• Determinant of the Jacobian: of the transformation - quantifies how much the image needed to be warped at that location to fit into the template space
• Voxel-based morphometry combines Jacobian and tissue segmentation images to compare gray matter volume in template space
  • Mechelli 2005 – Review of VBM
  • Antonopoulos 2023 – Recent comparison of VBM methods

What are fMRI data?

• fMRI measures brain activity over time by detecting changes in blood oxygenation.
• fMRI data are a time series of 3D brain images, collected every 1–3 seconds (TR = repetition time).
• Each voxel (3D pixel) contains a time series reflecting fluctuations in signal intensity across the scan.
• The signal includes BOLD (Blood Oxygenation Level-Dependent) response
  • Based on the fact that oxygenated and deoxygenated hemoglobin have different magnetic properties.
  • Neural activity, increased blood flow, local increase in oxygenated blood, measurable signal change.
• Data are collected while participants
  • Perform a task (task-based fMRI)
  • Rest with eyes open/closed (resting-state fMRI)
• fMRI is indirect: it measures blood flow, not neural activity itself.

Functional preprocessing pipeline

• Scientists have developed many software tools for preprocessing
• The analysis pipeline is very sophisticated
• Teams have developed pipeline tools that make it easier to reproducibly preprocess the data
  • Example CPAC functional pipeline

Functional pipeline

Slice Timing Correction

• fMRI data are collected sequentially in 2-d slices in the z-plane
  • Over the course of one TR (repetition time; approximately 2 seconds [0.5-3 seconds] ), the time the first slice and the last slice are acquired is almost the length of the TR
• Slice timing adjusts voxel intensities in each slice to a common acquisition reference (e.g., the first slice)
  
• Corrects timing offsets across the TR to align the hemodynamic response
  
• Uses interpolation (e.g., sinc, Fourier) based on slice order (alt‑z, seq, etc.)
• Multiband fMRI methods can acquire several slices simultaneously, reducing the TR.

Slice timing GIF in Andy’s Brain Book

Motion correction

• The participant will move across the scan time
• Motion correction in FSL
  1. Align volumes to mean image
• Motion correction in AFNI (another software)
  1. Align volumes to mean image
  2. Realign again to newly computed mean for improved accuracy

fsleyes data/RBC/HBN_BIDS/sub-NDARAA306NT2/ses-HBNsiteRU/func/sub-NDARAA306NT2_ses-HBNsiteRU_task-movieDM_bold.nii.gz

Distortion correction in FSL (skip this)

cd data/RBC/HBN_BIDS/sub-NDARAA306NT2/ses-HBNsiteRU/fmap
# merge the two phase encoding images
fslmerge -t AP_PA_merged.nii.gz sub-NDARAA306NT2_ses-HBNsiteRU_acq-fMRI_dir-AP_epi.nii.gz sub-NDARAA306NT2_ses-HBNsiteRU_acq-fMRI_dir-PA_epi.nii.gz
# write acquistion parameters x, y, z phase encoding directions. Final parameter is total readout time
echo "0 -1 0 0.04565\n0 1 0 0.04565" > acquisitionParams.txt
# apply topup
topup --imain=AP_PA_merged.nii.gz \
--datain=acquisitionParams.txt \
--config=b02b0.cnf \
--out=topup_AP_PA_b0 \
--iout=topup_AP_PA_b0_iout \
--fout=topup_AP_PA_b0_fout
# create mean of two unwarped images
fslmaths topup_AP_PA_b0_fout.nii.gz -abs -Tmean b0_mean.nii.gz
# brain extract the mean image
bet b0_mean.nii.gz b0_mean_brain.nii.gz -m
# convert Hz to radians
fslmaths topup_AP_PA_b0_iout.nii.gz -mul 6.283185307179586 field_rads.nii.gz

FEAT B0 options:

• Fieldmap image (rad/s) is field_rads.nii.gz
• Magnitude image is b0_mean_brain.nii.gz
• Effective EPI echo spacing (ms) - obtained from 4D fMRI data .json file "TotalReadoutTime": 0.0481411
• Unwarp direction - obtained from 4D fMRI data .json file "PhaseEncodingDirection": "j-",
• EPI TE not needed